clear;
clc;

%% Preset Parameters
% Benchmark Parameters
StdVec          =   [1 1 1 -1 -1 60 -10]'*0.01/4;
Param_Benchmark =   struct('KAPPA',0.00/4,'XI',1/82.5/4,'Pir_SS',0.02/4,...
                            'Flag_dE',0,'Flag_Price',0,...
                            'B_Total_SS',0.91,'b_lb',-0.29,...
                            'TrProb_H',0.96,'Prob_H',0.37,...
                            'TrProb_ext',0.92,'Prob_ext',0.33,...
                            'Fin_DomShare',0.79,'Fin_AdjCost',0.80,...
                            'Cons_FracH',0.60,...
                            'RHO_M_dom',0.68,'RHO_M_ext',0.81,'RHO_Y_H',0.50,'RHO_p_F',0.90,...
                            'Taylor_Pi',1.1,'Taylor_ir',0.87,...
                            'Cons_ElasTN',6.19,'Cons_ElasHF',6.19);
% International Integration Degrees
List_Prob_ext   =   (0.00:0.1:1.00)';
Num_Prob_ext    =   length(List_Prob_ext);

List_Fin_AdjCost=   [0,0.01,0.05,0.1,0.2,0.4,0.8,2,5,10,50,100,1000];
Num_Fin_AdjCost =   length(List_Fin_AdjCost);
%% Generate Different Results
SOLUTION        =   cell(Num_Prob_ext,Num_Fin_AdjCost);
for ii_Prob_ext=1:Num_Prob_ext
    TempParam       =   Param_Benchmark;
    TempParam.Prob_ext= List_Prob_ext(ii_Prob_ext);
    disp(['Solve the model :' num2str(ii_Prob_ext)]);
    % Solve the Steady State
    TempSolution    =   Solver_Main(Setup_PP(TempParam),StdVec);
    parfor ii_Fin_AdjCost=1:Num_Fin_AdjCost
        Temp_PP                 =   TempSolution.PP;
        Temp_PP.Fin_AdjCost     =   List_Fin_AdjCost(ii_Fin_AdjCost);
        Temp_EquJac             =   TempSolution.EquJac;
        Temp_EquJac.Portfolio   =   Equ_Portfolio(Temp_PP,TempSolution.SS,TempSolution.MODEL,1,[],[],1);
        Temp_fir_ord            =   JacAssemble(TempSolution.MODEL,Temp_EquJac);
        [gx,hx,gxhx_ExitFlag,gxhx_Info]=gx_hx(Temp_fir_ord,1);
        if gxhx_ExitFlag~=1
            warning('Perturbation Solution is problematic');
            continue;
        else
            Temp_gxhx           =   struct('gx',gx,'hx',hx);
        end
        Temp_IRF    =   IRF_1order(TempSolution.MODEL,Temp_gxhx,StdVec,[],4*20);
        
        SOLUTION{ii_Prob_ext,ii_Fin_AdjCost} ...
                    =   struct('PP',Temp_PP,'SS',TempSolution.SS,'MODEL',TempSolution.MODEL,'EquJac',Temp_EquJac,...
                               'fir_ord',Temp_fir_ord,'gxhx',Temp_gxhx,'IRF',Temp_IRF);
    end
end

SOLUTION_IntDegree          =   cell(Num_Prob_ext,1);
SOLUTION_IntDegree_AdjCost  =   cell(Num_Prob_ext,Num_Fin_AdjCost);

for ii_Prob_ext=1:Num_Prob_ext
    SOLUTION_IntDegree{ii_Prob_ext} ...
                    =   struct('PP',SOLUTION{ii_Prob_ext,1}.PP,...
                               'SS',SOLUTION{ii_Prob_ext,1}.SS,...
                               'MODEL',SOLUTION{ii_Prob_ext,1}.MODEL,...
                               'EquJac',SOLUTION{ii_Prob_ext,1}.EquJac);
    for ii_Fin_AdjCost=1:Num_Fin_AdjCost
        SOLUTION_IntDegree_AdjCost{ii_Prob_ext,ii_Fin_AdjCost} ...
                    =   struct('EquJac_Portfolio',SOLUTION{ii_Prob_ext,1}.EquJac.Portfolio,...
                               'IRF',SOLUTION{ii_Prob_ext,1}.IRF);
    end
end

save('Results_IntDegreeAdjCost.mat','SOLUTION','-v7.3');
% save('Results.mat','SOLUTION_IntDegree','SOLUTION_IntDegree_AdjCost','-v7.3');

%% IRFs
MatSize     =   size(SOLUTION);
IRFs        =   cell(size(SOLUTION));
for ii=1:prod(MatSize)
    [rr,cc]         =   ind2sub(MatSize,ii);
    IRFs{rr,cc}     =   SubFun_CollectIRF(SOLUTION{rr,cc});
end
save('IRFs.mat','IRFs','-v7.3');